
#include <libconfig.h>
#include <dirent.h>

#include "initialize_config.h"
#include "gsl/gsl_rng.h"
#include "gsl/gsl_randist.h"
#include "stdlib.h"
#include "math.h"
#include "global_definitions.h"
#include "string.h"

/*HELP FUNCTIONS DECLARATIONS*/
void fill_matrix_with_unique_random_values(gsl_matrix*,
                                           double (*pdf)(double, double),
                                           const gsl_rng*,
                                           const double,
                                           const int,
                                           const int);

int is_unique_element_in_column(const gsl_matrix* ,
                                const int ,
                                const double ,
                                const double );

int read_vector_from_cfg_file(const config_t*, const char* , gsl_vector**);



/*
 *INTERFACE FUNCTIONS
 */
int conf_read_simulation_config(const char* file_name, fsp_t* fsp)
{
    simp_t* sim_cfg = fsp->simp;
    FILE* cfg_file = NULL;
    config_t sim_config;
    int result = 0;
    char* message_buffer;

    memset(sim_cfg, 0, sizeof(simp_t));
    config_init( &sim_config );
    cfg_file = fopen(file_name, "r");
    if (cfg_file != NULL)
    {
        result = config_read( &sim_config, cfg_file );
        if (result == CONFIG_TRUE)
        {
            printf( "Confiuration file is succesfully read\n" );
            /*read the configuration*/
            result = config_lookup_int( &sim_config, PATH_NUM_OF_Z_ELECTRONS, &(sim_cfg->NumOfZElectrons) );
            result &= config_lookup_int( &sim_config, PATH_NUM_OF_TR_ELECTRONS, &(sim_cfg->NumOfTrElectrons) );
            result &= config_lookup_float( &sim_config, PATH_TEMPERATURE_X, &(sim_cfg->TemperatureX) );
            result &= config_lookup_float( &sim_config, PATH_TEMPERATURE_Y, &(sim_cfg->TemperatureY) );
            result &= config_lookup_float( &sim_config, PATH_TEMPERATURE_Z, &(sim_cfg->TemperatureZ) );
            result &= config_lookup_float( &sim_config, PATH_E0, &(sim_cfg->E0) );
            result &= config_lookup_float( &sim_config, PATH_B0, &(sim_cfg->B0) );
            result &= config_lookup_float( &sim_config, PATH_N0, &(sim_cfg->n0) );
            result &= config_lookup_float( &sim_config, PATH_PULSE_DURATION, &(sim_cfg->PulseDuration) );
            result &= config_lookup_float( &sim_config, PATH_PULSE_DELAY, &(sim_cfg->PulseDelay) );
            result &= config_lookup_float( &sim_config, PATH_PULSE_LENGTH, &(sim_cfg->PulseLength) );
            result &= config_lookup_float( &sim_config, PATH_PULSE_Z_SHIFT, &(sim_cfg->PulseZShift) );
            result &= config_lookup_float( &sim_config, PATH_DZ, &(sim_cfg->dz) );
            result &= read_vector_from_cfg_file( &sim_config, PATH_TIME_COLLECT_POINTS, &(sim_cfg->time_collect_points) );
            result &= read_vector_from_cfg_file( &sim_config, PATH_Z_COLLECT_POINTS, &(sim_cfg->z_collect_points) );
            result &= config_lookup_float( &sim_config, PATH_ABS_ERR, &(sim_cfg->abs_eps) );
            result &= config_lookup_float( &sim_config, PATH_REL_ERR, &(sim_cfg->rel_eps) );
            result &= config_lookup_string( &sim_config, PATH_SAVEDIR, &(sim_cfg->savedir) );
            result &= config_lookup_string( &sim_config, PATH_TEMPDIR, &(sim_cfg->tempdir) );
            result &= config_lookup_string( &sim_config, PATH_NAME_TEMPLATE, &(sim_cfg->name_template) );

        }
        else
        {
            message_buffer = config_error_file( &sim_config );
            printf("Error during passing configuration file:\n%s",
                   message_buffer);
            result = EXIT_FAILURE;
        }
        
    }
    else
    {
        printf("Configuration file could not be open\n");
        result = EXIT_FAILURE;
    }
    fclose(cfg_file);
    return result;
}

 int conf_generate_calculated_parameters(fsp_t* fsp,
            double (*time_envelope)(fsp_t*, double),
            double (*time_filling)(fsp_t*, double),
            double(*space_envelope)(fsp_t*, double))
 {
     fsp->csimp->w_ce = ELECTRON_CHARGE*(fsp->simp->B0)/(ELECTRON_MASS*LIGHT_SPEED);
     generate_initial_speed_distributions(fsp->simp, fsp->csimp);
     generate_pulse(fsp, time_envelope, time_filling, space_envelope);
 }


 /*
 *HELP FUNCTIONS DEFINITIONS
 */
 
int generate_initial_speed_distributions(simp_t* simp, csimp_t* csimp)
{
    int NumOfZSpeedEl = simp->NumOfZElectrons;
    int NumOfTrSpeedEl = simp->NumOfTrElectrons;
    const gsl_rng_type* T;
    gsl_rng* r;
    double sigmaX = 0;
    double sigmaY = 0;
    double sigmaZ = 0;
    double max_vz = 0;

    printf("Initial speed distribution generation is started\n");
    sigmaX = sqrt((((simp->TemperatureX)*1.602e-12)/ELECTRON_MASS));
    sigmaY = sqrt((((simp->TemperatureY)*1.602e-12)/ELECTRON_MASS));
    sigmaZ = sqrt((((simp->TemperatureZ)*1.602e-12)/ELECTRON_MASS));

    gsl_rng_env_setup();

    T = gsl_rng_default;
    r = gsl_rng_alloc(T);

    csimp->trsd = gsl_matrix_alloc(NumOfTrSpeedEl, 4);
    csimp->zsd = gsl_matrix_alloc(NumOfZSpeedEl, 2);

    fill_matrix_with_unique_random_values(csimp->zsd,  &gsl_ran_gaussian_pdf, r, sigmaZ, COLUMN_ZSPEED_PDF , COLUMN_ZSPEED_VALUE);
    fill_matrix_with_unique_random_values(csimp->trsd, &gsl_ran_gaussian_pdf, r, sigmaX, COLUMN_XSPEED_PDF , COLUMN_XSPEED_VALUE);
    fill_matrix_with_unique_random_values(csimp->trsd, &gsl_ran_gaussian_pdf, r, sigmaY, COLUMN_YSPEED_PDF , COLUMN_YSPEED_VALUE);

    max_vz = MAX(gsl_matrix_max(csimp->zsd), fabs(gsl_matrix_min(csimp->zsd)));
    csimp->dt = simp->dz/max_vz;
    
    gsl_rng_free(r);
    printf("Initial speed distribution generation is end\n");
    return EXIT_SUCCESS;
    
}

int generate_pulse(fsp_t* fsp,
        double(*time_envelope)( fsp_t*, double ),
        double(*time_filling)( fsp_t*, double ),
        double(*space_envelope)( fsp_t*, double )
)
{
    simp_t* simp = fsp->simp;
    csimp_t* csimp = fsp->csimp;
    double current_time = 0;
    double current_z = 0;
    double dt = csimp->dt;
    double dz = simp->dz;
    double pulse_duration = simp->PulseDuration;
    double pulse_delay = simp->PulseDelay;
    double pulse_length = simp->PulseLength;
    double pulse_shift = simp->PulseZShift;
    double amp = simp->E0;
    gsl_vector* pulse_time_envelope;
    gsl_vector* pulse_time_filling;
    gsl_vector* pulse_z_envelope;
    int num_of_t_points = ceil(simp->PulseDuration/csimp->dt);
    int num_of_z_points = ceil(simp->PulseLength/simp->dz);
    int i = 0;
    int j = 0;
    double p_0 = 0;
    double p_1 = 0;

    printf("Generating of pulse amplitude distribution is started\n");
    csimp->pulse = gsl_matrix_alloc( num_of_t_points, num_of_z_points );
    csimp->pulse_time = gsl_vector_alloc( num_of_t_points );
    csimp->pulse_z = gsl_vector_alloc( num_of_z_points );
    pulse_time_envelope = gsl_vector_alloc(num_of_t_points);
    pulse_time_filling = gsl_vector_alloc(num_of_t_points);
    pulse_z_envelope = gsl_vector_alloc(num_of_z_points);

    i = 0;
    for (current_time = pulse_delay; current_time < (pulse_delay + pulse_duration); current_time = i*dt + pulse_delay)
    {
        gsl_vector_set( csimp->pulse_time, i, current_time);
        gsl_vector_set( pulse_time_envelope, i, time_envelope(fsp, current_time));
        gsl_vector_set( pulse_time_filling,  i, time_filling(fsp, current_time));
        i++;
    }

    i = 0;
    for (current_z = pulse_shift; current_z < (pulse_shift + pulse_length); current_z = i*dz + pulse_shift)
    {
        gsl_vector_set( csimp->pulse_z, i, current_z);
        gsl_vector_set(pulse_z_envelope, i, space_envelope(fsp, current_z));
        i++;
    }
    
    for (i = 0; i < num_of_t_points; i++)
    {
        for (j = 0; j < num_of_z_points; j++)
        {
            gsl_matrix_set(csimp->pulse, i, j, amp*gsl_vector_get(pulse_time_envelope, i)*
                    gsl_vector_get(pulse_time_filling, i)*
                    gsl_vector_get(pulse_z_envelope, j));
        }
    }

    csimp->pulse_dz = gsl_matrix_alloc( num_of_t_points, num_of_z_points );
    gsl_matrix_set_zero( csimp->pulse_dz );
    for (i = 0; i < num_of_t_points; i++)
    {
        for(j = 1; j < num_of_z_points; j++)
        {
            p_0 = gsl_matrix_get(csimp->pulse, i, j - 1);
            p_1 = gsl_matrix_get(csimp->pulse, i, j);
            gsl_matrix_set( csimp->pulse_dz, i, j, E_M*(p_1 - p_0)/(dz));
        }
    }

    csimp->pulse_dt = gsl_matrix_alloc( num_of_t_points, num_of_z_points );
    gsl_matrix_set_zero( csimp->pulse_dt );
    for (i = 1; i < num_of_t_points; i++)
    {
        for(j = 0; j < num_of_z_points; j++)
        {
            p_0 = gsl_matrix_get(csimp->pulse, i - 1, j);
            p_1 = gsl_matrix_get(csimp->pulse, i, j);
            gsl_matrix_set( csimp->pulse_dt, i, j, E_M*(p_1 - p_0)/(dt));
        }
    }

    gsl_vector_free(pulse_time_envelope);
    gsl_vector_free(pulse_time_filling);
    gsl_vector_free(pulse_z_envelope);
    printf ("Generating of pulse amplitude distribution is succesfully end\n");
    return EXIT_SUCCESS;
}

int conf_free_configuration(fsp_t* fsp)
{
    simp_t* simp = fsp->simp;
    csimp_t* csimp = fsp->csimp;
    printf("Free configuration resources\n");
    gsl_vector_free(simp->time_collect_points);
    gsl_vector_free(simp->z_collect_points);
    gsl_matrix_free(csimp->zsd);
    gsl_matrix_free(csimp->trsd);
    gsl_matrix_free(csimp->pulse);
    gsl_matrix_free(csimp->pulse_dz);
    gsl_matrix_free(csimp->pulse_dt);
    gsl_vector_free(csimp->pulse_time);
    gsl_vector_free(csimp->pulse_z);
    printf("Configuration resources are free now\n");
    return EXIT_SUCCESS;
}

int read_vector_from_cfg_file(const config_t* config, const char* path, gsl_vector** out)
{
    int result = CONFIG_FALSE;
    config_setting_t* vector = NULL;
    vector = config_lookup(config, path);
    if (vector != NULL)
    {
        int idx = 0;
        /*determinate vector size*/
        for(idx = 0;  config_setting_get_elem(vector, idx); idx++);

        /*read values*/
        if (idx > 0)
        {
            int i = 0;
            (*out) =  gsl_vector_alloc(idx);
            for(i = 0; i < idx; i++)
            {
                gsl_vector_set(*out, i, config_setting_get_float_elem(vector, i));
            }
            result = CONFIG_TRUE;
        }
    }
    return result;
}

int is_unique_element_in_column(const gsl_matrix* matrix,
                                const int column,
                                const double element,
                                const double treshold)
{
    int i = 0;
    int rowsize = matrix->size1;
    int unique = 1;
    for (i = 0; (i < rowsize) && (unique = (fabs(element - (gsl_matrix_get(matrix, i, column)))) > treshold); i++);
    return unique;
}

void fill_matrix_with_unique_random_values(gsl_matrix* matrix,
                                           double (*pdf)(double, double),
                                           const gsl_rng* r,
                                           const double pdf_parameter,
                                           const int pdf_col,
                                           const int val_col)
{
    int NumOfRows = matrix->size1;
    int i = 0;
    double val_temp;
    while (i < NumOfRows)
    {
        val_temp = gsl_ran_gaussian(r, pdf_parameter);
        if (is_unique_element_in_column(matrix, val_col, val_temp, 0))
        {
            gsl_matrix_set(matrix, i, pdf_col, pdf(val_temp, pdf_parameter));
            gsl_matrix_set(matrix, i, val_col, val_temp);
            i++;
        }
    }
}